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Abstract. In this paper, we address the problem of local search for the 
falsification of hybrid automata with affine dynamics. Namely, if we are 
given a sequence of locations and a maximum simulation time, we re- 
turn the trajectory that comes the closest to the unsafe set. In order 
to solve this problem, we formulate it as a differentiable optimization 
problem which we solve using Sequential Quadratic Programming. The 
purpose of developing such a local search method is to combine it with 
high level stochastic optimization algorithms in order to falsify hybrid 
systems with complex discrete dynamics and high dimensional contin- 
uous spaces. Experimental results indicate that indeed the local search 
procedure improves upon the results of pure stochastic optimization al- 
gorithms. 
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1 Introduction 

Despite the recent advances in the computation of reachable sets in medium 
to large-sized linear systems (about 500 continuous variables) [UH], the verifi- 
cation of hybrid systems through the computation of the reachable state space 
remains a challenging problem jHlH]. To overcome this difficult problem, many 
researchers have looked into testing methodologies as an alternative. Testing 
methodologies can be coarsely divided into two categories: robust testing [5HZ] 
and systematic/randomized testing [SHTTj. 

Along the lines of randomized testing, we investigated the application of 
Monte Carlo techniques [T2] and metaheuristics to the temporal logic falsification 
problem of hybrid systems. In detail, utilizing the robustness of temporal logic 
specifications [T3] as a cost function, we managed to convert a decision problem, 
i.e., does there exist a trajectory that falsifies the system, into an optimization 
problem, i.e., what is the trajectory with the minimum robustness value? The 
resulting optimization problem is highly nonlinear and, in general, without any 
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obvious structure. When faced with such difficult optimization problems, one 
way to provide an answer is to utilize some stochastic optimization algorithm 
like Simulated Annealing. 

In our previous work |12) . we treated the model of the hybrid system as a 
black box since a global property, such as convexity of the cost function, cannot 
be obtained, in general. One question that is immediately raised is whether we 
can use "local" information from the model of the system in order to provide 
some guidance to the stochastic optimization algorithm. 

In this paper, we set the theoretical framework to provide local descent in- 
formation to the stochastic optimization algorithm. Here, by local we mean the 
convergence to a local optimal point. In detail, we consider the falsification prob- 
lem of afhne dynamical systems and hybrid automata with affinc dynamics where 
the uncertainty is in the initial conditions. In this case, the falsification problem 
reduces to an optimization problem where we are trying to find the trajectory 
that comes the closest to the unsafe set (in general, such a trajectory is not 
unique). A stochastic optimization algorithm for the falsification problem picks 
a point in the set of initial conditions, simulates the system for a bounded du- 
ration, computes the distance to the unsafe set and, then, decides on the next 
point in the set of initial conditions to try. Our goal in this paper is to provide 
assistance at exactly this last step. Namely, how do we pick the next point in 
the set of initial conditions? Note that we are essentially looking for a descent 
direction for the cost function in the set of initial conditions. 

Our main contribution, in this paper, is an algorithm that can propose such 
descent directions. Given a test trajectory s xo : R + i— > M. n starting from a point 
xq, the algorithm tries to find some vector d such that s Xo +d gets closer to 
the unsafe set than s Xo . We prove that it converges to a local minimum of the 
robustness function in the set of initial conditions, and demonstrate its advan- 
tages within a stochastic falsification algorithm. The results in this paper will 
enable local descent search for the satisfaction of arbitrary linear temporal logic 
specifications, not only safety specifications. 

2 Problem Formulation 

The results in this paper will focus on the model of hybrid automata with afhne 
dynamics. A hybrid automaton is a mathematical model that captures systems 
that exhibit both discrete and continuous dynamics. In brief, a hybrid automaton 
is a tuple 

H = (X, L, E, Inv, Flow, Guard, Re) 

where X C R" is the state space of the system, L is the set of control locations, 
E C L x L is the set of control switches, Inv : L — > 2 X assigns an invariant set to 
each location, Flow :LxI-> K™ defines the time derivative of the continuous 
part of the state, Guard : E — > 2 X is the guard condition that enables a control 
switch e and, finally, Re :Ix£->IxLisa reset map. Finally, we let 
H = L x X to denote the state space of the hybrid automaton H. 



Formally, the semantics of a hybrid automaton are given in terms of general- 
ized or timed transition systems |14j . For the purposes of this paper, we define 
a trajectory rjh starting from a point ho S H to be a function rjh ■ R+ — > H. 
In other words, the trajectory points to a pair of control location - continuous 
state vector for each point in time: % (t) = (l(t), s Xo (t)), where l(t) is the lo- 
cation at time t, and s xo (t) is the continuous state at time t. We will denote 
by loc(f7/j ) 6 L* U L u the sequence of control locations that the trajectory % 
visits (no repetitions). The sequence is finite when we consider a compact time 
interval [0, T] and r\ is not Zeno. 

Assumptions: In the following, we make a number of assumptions. First, 
we assume that for each location v £ L the system dynamics are affine, i.e., 
x — Flow(v, x) — Ax + b, where A and b are matrices of appropriate dimensions. 
Second, we assume that the guards in a location are non-overlapping and that 
the transitions are taken as soon as possible. Thirdly, we assume that the hybrid 
automaton is deterministic, i.e., starting from some initial state, there exists 
a unique trajectory rjh of the automaton. This will permit us to use directly 
results from [5]. We also make the assumption that the simulation algorithms 
for hybrid systems are well behaved. That is, we assume that the numerical 
simulation returns a trajectory that remains close to the actual trajectory on a 
compact time interval. To avoid a digression into unnecessary technicalities, we 
will assume that both the set of initial conditions and the unsafe set are included 
in a single (potentially different) control location. 

Let U C H be an unsafe set and let Djj '■ H \- > R + be the distance function 
to U, defined by 



where pr L is the projection to the set of locations, pv x is the projection to the 
continuous state-space and 



Definition 1 (Robustness) Given a compact time interval [0, T\, we define 
the robustness of a system trajectory ijh starting at some h = (l,x) G H to be 
f(h) = mmo< t <T Du(r]h(t)). When I is clear from the context, we'll write f(x). 

Our goal in this paper is to find operating conditions for the system which 
produce trajectories of minimal robustness, as they indicate potentially unsafe 
operation. This can be seen as a 2-stage problem: first, decide on a sequence of 
locations to be followed by the trajectory. Second, out of all trajectories following 
this sequence of locations, find the trajectory of minimal robustness. This paper 
addresses the second stage. The central step is the solution the following problem: 

Problem 1 Given a hybrid automaton %, a compact time interval [0,T], a set 
of initial conditions Hq C H and a point ho = (Iq,xo) €E Hq such that < 
f(ho) < +oo. find a vector dx such that h' — (Iq,xq + dx), loc(r)h ) = locirj^i ) 




du(x) = inf ||a; — u\ 



andf(h' )<f(h ). 



An efficient solution to Problem[T]may substantially increase the performance 
of the stochastic falsification algorithms by proposing search directions where the 
robustness decreases. In summary, our contributions are: 

— We formulate Problem U as a nonlinear optimization problem, which we 
prove to be differentiable w.r.t. the initial conditions. Thus it is solvable 
with standard optimizers. 

— We developed an algorithm, Algorithm[TJ to find local minima of the robust- 
ness function. 

— We demonstrate the use of Algorithm Q] in a higher-level stochastic falsifica- 
tion algorithm, and present experimental results to analyze its competitive- 
ness against existing methods. 

3 Finding a descent direction 

Consider an affine dynamical system in R™, 

x = F{x) = Ax + b 
which we assume has a unique solution 

s Xo (t) = e At x + c(t) 

where xq G Xq is the initial state of the trajectory 

Let U C K™ be the convex set of bad states, and U its closure. Note that even 
for linear systems, / : Ao M> M+ is not necessarily differentiable or convex. Our 
goal is to find the trajectory of minimum robustness. That is done by a local 
search over the set of initial conditions. 

Given an initial state Xq and a trajectory s Xo that starts at xo, define the 
time t* of its closest proximity to U, and the point u* G U which is closest to 
the trajectory: 

t* = arg min du(s Xo (t)) , u* — argmhi \\s Xo (t*) — u\\ 
t>o " ' ueu 

3.1 Partial descent based at the nearest point 

Given t* , choose an approach vector d! such that s XQ {t*) + dl is closer to U than 
Sx (t*)- Such a vector always exists given that s Xo has a positive distance to Li. 
Moreover, it is not unique. Thus we have 

f(x ) = \\s Xo (t*) - u*\\ > mhx\\s Xo (t*) + d' - u\\ 

U 

Define d = e~ At 'd'. Then 

f(x ) > min||s xo (t*) +d! - u\\ = min ||e A **a;o + c(t) + e At *d - u\\ 

U U 

> min min || (x + d)e At + c(t) - u\ \ = f(x + d) > 



and d is a descent direction, provided that xq + d G Xo. 
It is easy to see that for any xo G Xo and d' G R", 

Sx +e-^*d'(*) = s xa {t*)+d' 

so the new distance is achieved at the same time t* as the old one. This new 
distance du(s Xo +d{t*)) is an upper bound on the new trajectory's robustness. In 
general, the new trajectory's robustness might be even smaller, and achieved at 
some other time t' ^ t* . 

As pointed out earlier, the approach d' is not unique. The requirement on 
d' is that s Xo (t*) + d' be closer to U than s XQ (t*). So define the set P(xo]t*) of 
points that are closer to U than s Xo (t*) (see Fig. [1]): 

P(x Q ;t*)^{xeR n \du(x)<f(x )} (1) 



f P(x ;t*) N 
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Fig. 1. The unsafe set U and the set P(xo;t*). The system trajectory s xo appears as 
a dashed curve. 

Then d! must satisfy s Xo {t*)+d' G P(x ;t*) ode e- At '(P(x Q ; t*)-s Xo (t*)). 
Combined with the requirement that xq + d G Xq , we get 

d G (Xq — x ) p| e~ At * [P(x ; t*) - s Xo (t*)} 

Any point in the above descent set is a feasible descent direction. As a special 
case, it is easy to verify that d = e~ At (u — s Xo (t*)), for any u G U, is a descent 
direction that leads to robustness. Coupled with the requirement that xo + d 
must be in Xo, it comes 

de (X Q ~x )f]e- At '(U~s Xo (t*)) 

If computing P is too hard, we can approximate it with the following IA U : 
imagine translating U along the direction v = s Xo (t*) — u* , so it is being drawn 
closer to s Xo (t*), until it meets it. Then we claim that the union of all these 
translates forms a set of points closer to U than s Xo (t* ) : 



Proposition 1. Let IA be a convex set, s Xo (t*) a point outside it, and IA V (v) 
be the Minkowski sum of IA and {av\a £ [0,1]}. Then for any p in v). 
duip) < d u (s X0 (t*)) 

Proof. U u is convex by the properties of Minkowski sums. Let u £ dU. Then for 
any a < 1, djjiu-\- av) < \\u + av — u\\ — a\\v\\ = af(xo) < J(xq). So translates 
of boundary points are closer to U than s xo (t*). 

Now we show that all points in IA U /U are translates of boundary points. 
Consider any point p = u + av in IA U /IA: u is in IA, but p is not, so the line [u,p] 
crosses dlA for some value a°: u + a°v £ dlA. And, p = u + av = (u + a°v) + 
(a — a°)v, so by what preceded, du(p) < f{x). 
When p £ IA, of course, du{p) = < f(x). 

We have thus defined 3 possible descent sets: U C IA U C P{xq; t*). 
3.2 Implementation 

The question we address here is: how do we obtain, computationally, points in 
the descent set W, where W = U,P(xq) or U u (v)l The following discussion is 
based on Chapters 8 and 11 of [To] . 

Since we're assuming Xq and U to be convex, then the descent set is also 
convex. Describe Xq with a set of Nx inequalities qiix) < where the qi are 
convex and differentiable, and W = {x\pi{x; xq) < 0,i = l...k} for convex dif- 
fcrcntiable pi (the particular form of the pi will depend on the descent set at 
hand). We assume dom pi — dom qi =R n . 

Given an already simulated trajectory s Xo and its time of minimum robust- 
ness t* , we are looking for a feasible x\ such that s Xl (t) £ W for some t. Thus 
we want to solve the following feasibility problem 

min v 

s.t. pi{s x {t);x Q )<u,i = l...k (t-PDP(xo)) ( 2 ) 
qi(x) <v,i = l...N x 

This is a convex program, which can be solved by a Phase I Interior Point 
method [15] . A non-positive minimum v* means we found a feasible x; if W = IA, 
then our work is done: we have found an unsafe point. Else, we can't just stop 
upon finding a non-positive minimum: we have merely found a new point X\ 
whose robustness is less than xo's, but not (necessarily) 0. So we iterate: solve 
t-PDP(io) to get x\, solve t-PDP(ii) to get x 2 , and so on, until f(xi) — 0, a 
maximum number of iterations is reached, or the problem is unsolvable. If the 
minimum is positive, this means that for this value oft, it is not possible for any 
trajectory to enter IA at time t. 

The program suffers from an arbitrary choice of t. One approach is to sample 
the trajectory at a fixed number of times, and solve ^ for each. This is used in 
the experiments of this section. A second approach, used in the next section, is 
to let the optimization itself choose the time, by adding it to the optimization 
variable. The resulting program is no longer necessarily convex. 



3.3 Numerical Experiments 



In this section, we present some numerical experiments demonstrating the prac- 
tical significance of the previous theoretical results. 

Example 1 We consider the verification problem of a transmission line U 6] /. 
The goal is to check that the transient behavior of a long transmission line has 
acceptable overshoot for a wide range of initial conditions. Figure fJI shows a 
model of the transmission line, which consists of a number of RLC components 
(R: resistor, L: inductor and C: capacitor) modeling segments of the line. The left 
side is the sending end and the right side is the receiving end of the transmission 
line. 




Fig. 2. RLC model of a transmission line. 



The dynamics of the system are given by a linear dynamical system 

x{t) = Ax(t) + bV ln (t) and V out (t) = Cx(t) 

where x(t) G M 81 is the state vector containing the voltage of the capacitors and 
the current of the inductors and Vi n (t) £ M is the voltage at the sending end. 
The output of the system is the voltage V ou t (t) £ K at the receiving end. Here, 
A, b and C are matrices of appropriate dimensions. Initially, we assume that 
the system might be in any operating condition such that x(0) £ [— 0.1,0. 1] 41 x 
[— 0.01, 0.01] 40 . Then, at time t = the input is set to the value Vi n (t) = 1. 

The descent algorithm is applied to the test trajectory that starts from x(0) = 
and it successfully returns a trajectory that falsifies the system (see Fig. [#[). 

4 Hybrid systems with affine dynamics 

We now turn to the case of hybrid systems with affine dynamics in each location. 
The objective is still to find a descent direction in Ho, given a simulated trajec- 
tory rjh originating at point h$ £ Hq. Note that since we have assumed that 
pr L (_ffo) is a singleton set, the problem reduces to finding a descent direction in 

X Q = P r x(#o)- 

Assumptions. At this point, we make the following assumptions: 
a. The continuous dynamics in each location are stable^ 

1 This is not a restrictive assumption since we can also consider incrementally stable 
systems [17j . and even unstable linear systems [18] . 
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Fig. 3. The unsafe set U , the initial test trajectory starting from x(0) = and the 
trajectory that falsifies the system. 



b. For every transition e G L 2 , the resets Re(-, e) are diffcrcntiable functions 
of their first argument. 

c. Conditions 4 and 5 of Theorem III. 2 in [19] are satisfied, namely: for 
all i, there exists a differentiable function Oi : K" h-> K such that Inv{li) = 
{x G R n |<7i(a;) > 0}; and, for all z,x such that (Ti(x) = 0, the Lie derivative 
Lpai(x) =fi 0. This allows us to have a differentiable transition time t x of the 
trajectory starting at the initial point x £ Xq. 

d. The sequence of locations loc(r)} l0 ) enters the location of the unsafe set. 
This is required for our problem to be well-defined (specifically, for the objective 
function to have finite values). The task of finding such an h is delegated to the 
higher-level stochastic search algorithm, within which our method is integrated. 

4.1 Descent in the Robustness Ellipsoid 

Consider a trajectory r)h with positive robustness, with loc(?7fc, ) = Zo^i • • ■ In- 
This is provided by the simulation. Let the initial set Xq be in location Iq and 
let lu denote the location of U. In order to solve Problem [TJ we assume that 
lu appears in loc(?7fc, ) (see Assumption d above) - otherwise, f(ho) = +oo and 
the problem as posed here is ill-defined. We search for an initial point h' G Ho 
(actually x' G Xq), whose trajectory gets closer to the unsafe set than the 
current trajectory rjh ■ 

In order to satisfy the constraints of Problem [lj we need to make sure that 
the new point h' that we propose generates a trajectory that follows the same 
sequence of locations as rjh ■ This constraint can be satisfied using the notion 
of robust neighborhoods introduced in [B]. In [BJ, it is shown that for stable 
systems and for a given safe initial point hg = (Iq,xq), there exists an 'ellipsoid 
of robustness' centered on xq, such that any trajectory starting in the ellipsoid, 




remains in a tube around T]h . The tube has the property that all trajectories in 
it follow the same sequence of locations as % . Therefore, we restrict the choice 
of initial point to Xq f] E(xo), where E(y) — {x\(x — y) T R^ 1 (x — y) < 1} is the 
ellipsoid of robustness centered on xq, with shape matrix R. Formally, in [B], the 
following result was proven. 

Theorem 1 Consider a hybrid automaton %, a compact time interval [0,T], 
a set of initial conditions Hq C H and a point ho = (lo,xo) G Hq. Then, we 
can compute a number e > and a bisimulation junction <f){x\,X2) = {x\ — 
X2) 7 M {x\ — X2), where M is a positive semidefinite matrix, such that for any 
x' e {y 6 X (j>(x ,y) < e}, we have loc(n ho ) = loc(n^ ^ o) ). 

Remark 1 (i) In Jfjjj, in the computation of e, we also make sure that any point 
in the robust neighborhood generates a trajectory that does not enter the unsafe 
set. In this work, we relax this restriction since our goal is to find a point that 
generates a trajectory that might enter the unsafe set. (ii) In view of Theorem 
[7J the shape matrix for the ellipsoid is defined as R — e 2 M~ 1 . 

We now proceed to pose our search problem as a feasibility problem. Let to 
be the time at which s Xo is closest to U. We choose P(xq; to) as our descent set: 
recall that it is the set of all points which are closer to U than s X(i (to) (Def. [T]). 
Therefore, if we can find x* £ Xq C\ E(xq) such that s x * (t*) £ P(xq; to) for some 
t*, it follows that f(x*) < f(x ). To simplify notation, let W = P(x ;t ) be the 
descent set. As before, it is assumed that W = {x <E R n : Pi(x) < 0,i = 1 . . . k} 
for differentiable pi. The search problem becomes: 

Given n ho , find x* £ X f) E(x ) and t* > 0, such that s x *{t*) E W. This is 
cast as an optimization problem over z e R" x R + x M: 

min v 

z=(x,t,i/) 

s.t. Cqx — go < 

(x — Xq) T P^ 1 (x — Xq) — 1 < V 

Pi(s x (t);x ) <v,i = l...k 
where s x (t) = W x (V(i ,x)(t)) and X = {x\C Q x ~ go < 0}. 

Remark 2 Note that Problem Q is specific to a choice of initial point xo; this 
will be important in what follows. In our implementation, the first constraint is 
specified as bounds to the optimization and so is always satisfied. 

Later in this section, we discuss how to solve this optimization problem. For 
now, we show how solving this problem produces a descent direction for the 
robustness function. For convenience, for z = (x,t,v), we define the constraint 



(3) 



functions 



Gq(z) = Cox ~ g (4a) 
G E (z) = {x-x ) T p-\x-x Q )-l (4b) 
/ 'pi(s x (t);x Q y 

(4c) 



Gw(z) 



A point z is feasible if it satisfies the constraints in Problem ([3]). Finally, define 
the objective function F(z) = v. 

The objective function F(z) measures the slack in satisfying the constraints: 
a negative v means all constraints are strictly satisfied, and in particular, Gyy. 
Thus, we have a trajectory that enters W and, hence, gets strictly closer to U. 
This reasoning is formalized in the following proposition: 

Proposition 2. Let z* = (x* ,t* ,v*) be a minimum of F(z) in program 
Then f(l ,x*) < /(Z ,x ). 

Proof. It is assumed that the optimizer is iterative and that it returns a solution 
that decreases the objective function. In what follows, for a vector y G R™, maxy 
is the largest entry in y. 

We first remark that for a given x and t that satisfy the constraints in ([3]), 

z = (x, t, maxjG^a;, t), Gw(x, t)}) 

is feasible, and F(z) < F(x, t, v) for any feasible (x, t, v). Therefore, we may only 
consider points with F(z) — v = max{G£;(x, t), Gw(x, t)}. 

Let zq = (xq, to,i>o) be the initial point of the optimization. Because xq is the 
center of E(xq), Ge(zo) — — 1. And, because s Xo (t ) € dW, maxGw(zo) = 0. 
Thus vq — 0. Therefore, at the minimum z* — (x*,t*,v*) returned by the 
optimizer, v* < vq — 0. In particular, Gyv(z*) < 0, and the new trajectory s x * 
enters W. Therefore, its robustness is no larger than that of the initial trajectory 
s Xo - 

We now address how Problem |3] might be solved. Functions F, Gq and 
Ge are differentiable in z = (x,t,is). It is not clear that Gyv, or equivalently, 
Pi(s x (t); xq), as a function of z, is differentiable. We now show that under some 
asumptions on the pi, for trajectories of linear systems, pi is in fact differentiable 
in both x and t, over an appropriate range of t. This implies differentiability in 
z. Therefore, standard gradient-based optimizers can be used to solve Problem 

m 

For the remainder of this section, we will re- write s x (t) as s(x, t) to emphasize 
the dependence on the initial point x. s^(x, r) will denote the point, at time r, 
on the trajectory starting at x G Inv(li), and evolving according to the dynamics 
of location i. When appearing inside location-specific trajectories such as sW ; the 
time variable will be denoted by the greek letter r to indicate relative time: that 



is, time measured from the moment the trajectory entered U, not from datum 
0. s(x,t) (without superscript) will denote the hybrid trajectory, traversing one 
or more locations. We will also drop the £0 from pi(y;xo), and write it simply 
as Pi(y). 

We first prove differentiability in x. Therefore, unless explicitly stated oth- 
erwise, the term 'diffcrcntiable' will mean 'diffcrcntiable in x\ Start by noting 
that p i (s{x,t)) is a composite function of x £ X . Since pi is differentiable, it 
is sufficient to prove that s(x,t) is differentiable. The hybrid trajectory s(x, •) is 
itself the result of composing the dynamics from the visited locations Iq, Zjv-i- 

Recall that E(xq) is the ellipsoid of robustness centered at xq. As shown by 
Julius et al. [6], the following times are well-defined: 

Definition 2 (Transition times) Given xq S Xq, let Eq = 'mt(E(xo) |~) Xq). 
ti is the time at which trajectory s(xq) transitions from Inv(li—\) into Inv(li) 
through guard Guard{li-i,U). 

is the maximal time for which the image of Eq under the hybrid dynamics is 
contained in Inv(li-\): 

tj = max{i|s(£'o, t) C 

In other words, t~ is the time at which occurs the first U-i-to-U transition of a 
point in s(Eq). 

tf is the minimal time for which the image of Eq under the hybrid dynamics 
is contained in Inv(k): 

t+ = mm{t\s(E ,t) C Inv(k)} 

In other words, tf is the time at which occurs the last li-\-to-li transition of a 
point in s(Eq). 

For a given point x £ Xq, t % ~ 1 ^' 1 (t^T 1 ^' 1 ) is the absolute (relative) transition 
time of trajectory s^' (x) from Inv{U-i) into Inv{h) through guard Guar d(h-i, h) . 
Thus, for example, t\ — t^ 1 — T®^* 1 and t 2 = t 1 ^ 2 = t®^ 1 + t}^ 2 , with 
2/o = s°(xq, t®^ 1 ). When the transition is clear from context, we will simply 
write t x (t x ). 

We will first show differentiability of a trajectory that visits only 2 locations 
Zo and 1%: 

s(x ,t) = 3 W(Re(sW(x ,t x ),(l ,h)),t -tx) (5) 

Example 2 We first present a simple ID example to illustrate the definitions 
and the idea of the proof. Consider the hybrid system with three locations 

U = (M, {0, f , 2}, {(0, f ), (f , 2)}, Inv, Flow, Guard, Id) 

where Inv(l) = M for 1 = 0,1, 2, and the flow is defined by 



Flow(l, x) = x(t) 



x(t) if I G {0, 2} 
-x(t) if 1 = 1 



The guards are Guard(0,l) = {1} and Guard(l,2) — {1/4}. Id is the identity 
map, so there are no resets. The initial set is Xq = [0, 1/2]. The solutions in the 
individual locations are then 

s m (x,t)=e t x 

s (1) 0,i) = e^x 
8 W(x,t) = e l x 

We can solve, in this simple case, for t^ 1 : e Tm x = 1 =>■ t x ^ 1 = ln(l/x). 
Similarly for t x ^ 2 : e~ T * • 1 = 1/4 t x ^ 2 = ln(Ax). 

We first show differentiability of the trajectory over locations and 1. We 
then do the same for a trajectory over locations 1 and 2. Then we stitch the 
two together and show differentiability over 3 locations. For locations and 1: 
s(x,t) = s^(s^(x,t x ),t - t x ) = s«(l,f - t x ) = e -(*-**) • 1 = e -*/x 
Jjs(x,t) = — ^r- 

Moving on the trajectory over locations 1 and 2, the procedure is the same: 
from an initial point x G Guard(0,l) = {1}, for a fixed (relative time) t G 
(t 2 -t u t- -h): s{x,t) = S M{sW(x,T x ),T-T x ) = S ( 2 )(l/4,r + ln(l/4x)) = 

Finally we stitch up the 2 portions of the trajectory: x G Xq, t G [^j^]- 
s(x,t) = s( 2 H s ( 1 H s ( a \x,t 1 ),t 2 - tl ),t-t 2 ) = aCV^Ma-tiJ.t-ta) = 
s( 2 )(l/4,t - t 2 ) = e^/A. Since t 2 = t ^ 1 + r\^ 2 = ln(l/a;) + ln(4 • 1) = 
In(4/a;) s(x,t) = ^ e ^ x /^ = xe t /16 f t s{x,t) = e*/16. 

We now prove the general case. 

Proposition 3. Let xq G Eq, and fix t G {ti,t 2 \. Consider the hybrid trajec- 
tory over 2 locations in Eq.§5$. If Assumptions a-d are satisfied, then s(x,t) is 
differ entiable at Xq . 

Proof. In what follows, e = (Iq,Ii)- 

s^(x,t x ) = e T * A °x + I e (T *- s)Ao bds 







e x + e 




terml term2 

term3 



Terms 1 and 2 are clearly differentiable in x. For term3, write M(t) = L e sA °bds 
so term3 = M(t x ). M(t) is differentiable by the 2 nd Fundamental Theorem of 
Calculus and its derivative is M'(t) = e~ tA °b. As a consequence of Assump- 
tion c, t x itself is differentiable in x (Lemma III. 3 in |19|). and the chain 
rule allows us to conclude that term3 is differentiable in x. Thus s^(x,t x ) 
is differentiable over Eq. Since Re(-,e) is differentiable by Assumption b, then 



i£e(s(°) (x, t x ), e) is differentiable over £?o- Note that Eq is open and s^ ^ is con- 
tinuous, so U = {w G W n \w = s^°'(x,t x ) for some x G Eq} C Guard(e) is open. 
Since Re{-,e) is continuous, then Re(U,e) is open. Next, 



Using the same argument as above, terms 4, 5 and 6 are differentiable in x. In 
conclusion, s{x, t) is differentiable at over Eq, and this ends the proof. 

The following proposition generalizes Prop. [3] to trajectories over more than 
2 locations. 

Proposition 4. Fix t € (tjv"_i,T], and consider the hybrid trajectory over N > 
1 locations. Then s(x,t) is differentiable at xq for all xq G Eq. 

Proof. We argue by induction over the number of locations N. The base case 
N — 1 is true by hypothesis, and the case N — 2 has been proven in Prop. 
[3] For N > 2 and £ < ijv-i, let CC^i*) be the trajectory over the first ./V — 1 
locations, so that s(x, t) = s'- N ~ 1 \Re(((x, rjy-a)) {In -2, In-i)), t— By the 
induction hypothesis, C{ x ^) 1S differentiable at xo- Then £ and s^" 1 ) satisfy 
the conditions of the case N = 2. 

Differentiability with respect to time is easily proven: 

Proposition 5. Let xq 6 Eq and t G {tjsr—i,T), that is, a time at which the 
trajectory is in the last location. Consider the hybrid trajectory over N > 1 
locations. Then s{xo,t) is differentiable in t over [£jv_i,T). 

Proof. s{xo,t) = s( Ar_1 )(.x , t — tjv-i)' The location-specific trajectories sW(a;, •) 
are solutions of differential equations involving at least the first time derivative. 
Therefore, they are smooth over {tN-i,T). This implies differentibility of the 
hybrid trajectory s{xq, •) over the same interval. At t = T, the trajectory is only 
left-differentiable, since it's undefined from the right. 

The following result is now a trivial application of the chain rule to pi o s: 

Proposition 6. Let xq G Eq, t G (t/v-i,? 1 ). If Pi is differentiable for all i = 
1, . . . ,k, then Gyv *s differentiable in z over Eq x R + x R. 

We choose Sequential Quadratic Programming (SQP), as a good general- 
purpose optimizer to solve Problem [3] SQP is a Q-quadratically convergent it- 
erative algorithm. At each iterate, Gw{xi 1 ti 1 Vi) is computed by simulating the 
system at a;,. This is the main computational bottleneck of this method, and 
will be discussed in more detail in the Experiments section. 



s{x, t) = s« {Re{s^{x, t x ), e), t - t x ) 




term6 



Algorithm 1 Robustness Ellipsoid Descent (RED) 



Input: An initial point xo £ Xo, and corresponding to- 
Output: zq. 
1: Initialization: i = 

2: Compute z* = (x* ,t* = minimum of Prot|3lWi]. 
3: while < do 

4: «- x* 

5: = argmint du(sx i+1 {t)) 

6: W 1+ i = P(x i+ i) 

7: Compute z* = (x* ,t* ,v*) — min of ProtJ3JWi+i]. 
8: i = i + l 
9: end while 
10: 

11: Return zq = z* 



4.2 Convergence to a local minimum 

Solving Problem §),fora given W, produces a descent direction for the robust- 
ness function. However, one can produce examples where a local minimum of 
F(-) is not a local minimum of the robustness function /. This section derives 
conditions under which repeated solution of Problem yields a local minimum 
of the robustness function. 

For i = 0, 1, 2, . . . , let Xi £ X n f] E(xi-i), and let U be the time when s Xi 
is closest to U. Let Wj = P(xi\ti) be the descent set for this trajectory. For 
each Wj, one can setup the optimization Problem ([3]) with W = Wi, and initial 
point (xi,ti,0); this problem is denoted by Prot|3[Wi]. (Recall from the proof 
of Proposition [5] that v — at the initial point of the optimization problem). 
Finally, let z* = (x*,t*,v*) be the minimum obtained by solving Prot|3[ W<]. 

Algorithm[T]describes how to setup a sequence of optimization problems that 
leads to a local minimum of /. It is called Robustness Ellipsoid Descent, or RED 
for short. 

Proposition 7. Algorithm^ (RED) terminates 

Proof. Proposition[2]holds for each problem ProtOJWi]. Therefore, each solution 
with Vi < gives a trajectory s x * with a smaller robustness than s x * i : f{x*) < 
f( x i-i)- Thus (/(xj))jgjv is a decreasing sequence, lower bounded by 0. There- 
fore, it converges to a limit r > 0. But how to prove that this limit is indeed a 
minimum of f ? 

Proposition 8. Assume that Algorithm [7] halts at a point zq = (xQ,tQ,VQ), 
for which there exist t\ , t% such that: 

- < ti < t Q < h 

- d u (s x _ Q {t)) > f(x Q )Vt eT R ^ [0,ii] U [t 2 ,T), and 

- <2 — t\ is 'sufficiently small'. 



Then xq is a local minimum of the robustness function f . 



Proof. We assume that the trajectory starting at xq is safe - otherwise, we're 
done since we found an unsafe trajectory. 

Two tubes will be constructed: one contains s XQ over (iijia)) the other con- 
tains it over Tr. They are such that no trajectory in them gets closer to Wq 
than s XQ . Then it is shown that all trajectories in a neighborhood of xq are 
contained in these tubes, making xq a local minimum of the robustness function 
/■ 




Fig. 4. [Proof of Prop[8] All trajectories starting in a neighborhood of xq will be 
contained in the orange tube over Tr and in the green tube over (ti,t2) 



By the halting condition, vq — 0. Since the optimizer always returns a local 
minimum of the objective function F, there exists a neighborhood N(zq) of zq 
such that for all z e N(z Q ),F(z) > F(z Q ) =v Q =Q 

V(x, t, v) E N(z Q ), s x (t) £ intWq 

o V(as, t, v) e N(zq), du(s x (t)) > f{x Q ) 
N(zq) can be expressed as 

n ( z q) = B(x Q ,e) x (t 3 ,t 4 ) x {-v u vi) 

e > 0, vi > 0, B(x Q , e) C E(x ) n X 

(Since M™ x M + x 1 is a finite product, the box and product topologies are 
equivalent, so it doesn't matter which one we use.) 

We now precise the notion of 'small enough': we require that 

(iiM C (ts,U) 



Therefore 



Mx G B(x Q ,e),te (ti,t 2 ),d u (s x (t)) > f(x Q ) 

Thus 

WxeB(x Q ,e), inf d u {s x (t)) > f{x Q ) (6) 

te(ti,ta) 

We now study the behavior of trajectories starting in B(xQ,e) over the re- 
maining time periord Tr. Recall that £Y C Wq. Let w° be any point on the 
boundary <9Wq. Then 

Vi G T R ,d u {s XQ (t)) > f(x Q ) = d u {w°) > 
=^Vi GTR,d WQ (s XQ (t)) >0 

Then 

yl = inf{d WQ ( S;EQ (t))|teT fl }>0 

s x is continuous as a function of x for every t, therefore 

36 > s.t. a: G B(x Q ,S) d(s XQ (t),s x (t)) < A 

Pick any point tu G Wq. Then Vx G B(xq,6) and t E T R 

d(s XQ (t),w) < d(s XQ (t),s x (t)) +d(s x (t),w) 
<A + d(s x (t),w) 
=>■ d(s x (t),w) > d(s XQ (t), to) - A 

Minimizing both sides over w G Wq, 

d WQ (s x (t)) > dw Q {s XQ (t)) - A > 
^ inf d Wo (»x(t)) >0 

In conclusion 

Vx G B(iq,<5), inf > f(x Q ) (7) 

Putting Eqs.(H]) and © together, it comes that Vx G B{xq, min{e, 8}) 

inf d u (s x (t)) > f(x Q ) 

t£K_l_ 

^V.t G fl(afQ,min{c,*}),/(x) > f(x Q ) 
and xq is a local minimum of the robustness /. 



Algorithm 2 RED with Simulated Annealing (SA+RED) 
Input: An initial point x G Xq. 
Output: Samples & C Xq. 

Initialization: BestSoFar = x, /{, = /(BestSoFar) 
1: while f(x) > do 
2: x 1 — ProposalScheme(a;) 
3: a = exp (-f3(f(x') - /„)) 
4: if [7(0, 1) < a then 
5: a;* = RED(as') 

6: x = x* 

7: else// Use the usual acceptance criterion 

8: a = exp (-/3(/(x') - /(x))) 

9: if U(0, 1) < a then a; = x' 

10: end if 

11: end if 

12: (BestSoFarJf,) = BetterOf(a;, BestSoFar) 
13: end while 



4.3 Ellipsoid Descent with Stochastic Falsification 

As outlined in the introduction, the proposed method can be used as a sub- 
routine in a higher-level stochastic search falsification algorithm. A stochastic 
search will have a ProposalScheme routine: given a point x in the search space, 
ProposalScheme will propose a new point X clS cL falsification candidate. Robust- 
ness Ellipsoid Descent (RED) may then be used to further descend from some 
judiciously chosen proposals. Algorithm [2] illustrates the use of RED within the 
Simulated Annealing (SA) stochastic falsification algorithm of [12]. {7(0, 1) de- 
notes a number drawn uniformly at random over (0, 1). Given two samples x and 
y, BetterOf(a;, y) returns the sample with smaller robustness, and its robustness. 

For each proposed sample x', it is attempted with certainty if its robustness 
is less than the smallest robustness fa found so far. Else, it is attempted with 
probability e -^ f ^- f ^ (lines 3-4). If x' is attempted, RED is run with x' as 
starting point, and the found local minimum is used as final accepted sample 
(line 6). If the proposed sample is not attempted, then the usual acceptance- 
rejection criterion is used: accept x' with probability min{l, e~^^ x )~/( a: ))}. As 
in the original SA method, ProposalScheme is implemented as a Hit-and-Run 
sampler (other choices can be made). The next section presents experimental 
results on three benchmarks. 

4.4 Experiments 

This section describes the experiments used to test the efficiency and effective- 
ness of the proposed algorithm SA+RED, and the methods compared against 
it. 

We chose 3 navigation benchmarks from the literature: NavO (4-dimensional 
with 16 locations) is a slightly modified benchmark of [20], and it is unknown 



Fig. 5. The navigation benchmark example. 



whether it is falsifiable or not. Navl and Nav2 (4-dimensional with 3 locations) 
are the two hybrid systems in the HSolver library of benchmarks [H) , and are 
falsifiable. We also chose a filtered oscillator, Fosc (32-dimensional with 4 lo- 
cations), from the SpaceEx library of benchmarks [22]. We describe the NavO 
benchmark that we used, as it a slightly modified version of the benchmark 
in [20]. 

Example 3 (Navigation Benchmark [20]) The benchmark studies a hybrid 
automaton H. with a variable number of discrete locations and 4 continuous vari- 
ables x\, X2, y%, 2/2 that form the state vector x — [x% X2 yi y2\ T ■ The structure 
of the hybrid automaton can be better visualized in Fig. [5] The invariant set of 
every location is an 1 x 1 box that constraints the position of the system, 

while the velocity can flow unconstrained. The guards in each location are the 
edges and the vertices that are common among the neighboring locations. 

Each location has affine constant dynamics with drift. In detail, in each lo- 
cation (i,j) of the hybrid automaton, the system evolves under the differential 
equation x ~ Ax — Bu(i,j) where the matrices A and B are 



A = 



1 
1 

-1.2 0.1 

0.1 -1.2 



id B = 





-1.2 0.1 
0.1 -1.2 



and the input in each location is 

u(i,j) = [sin(7rC(i, j)/4) cos(irC(i, j)/4)] T . 

The array C is one of the two parameters of the hybrid automaton that the user 
can control and it defines the input vector in each discrete location. Here, we 
consider the input array denoted in Fig. [5] 



The set of initial conditions is the set Hq = {13} x [0.2 0.8] x [3.2 3.8] x 
[—0.4 0.4] 2 (green box in Fig. [3J) and the unsafe set isU = {4} x {i £ I 4 | ||x — 
(3.5 0.5 0)|| < 0.3} (red circle in Fig. [5|). This is slightly modified from the 
original benchmark to simplify the programming of the pi functions. Sample tra- 
jectories of the system appear in[5\for initial conditions [0.8 3.2 — 0.2 0.35] T 
(red trajectory) and [0.4 3.3 — 0.1 — 0.1] T (blue trajectory). Note that the two 
trajectories follow different discrete locations. 

The methods compared are: SA+RED, pure Simulated Annealing (SA) [12] . 
mixed mode-HSolver (mm-HSolver) |21j , and the reachability analysis tool SpaceEx 
[22 . Because ours is a falsification framework, SpaceEx is used as follows: for 
a given bound j on the number of discrete jumps, SpaceEx computes an over- 
approximation R(j) of the set reachable in j jumps R(j) : R{j) C R(j). If 
R(j) C\U is empty, then a fortiori R(j) C\U is empty, and the system is safe if 
trajectories are restricted to j jumps. If, however, R(j) flU + 1 0, no conclusion 
can be drawn. 

Because SA and SA+RED are stochastic methods, their behavior will be 
studied by analyzing a number of runs. A regression will mean a fixed number of 
runs, all executed with the same set of parameters, on the same benchmark, mm- 
HSolvcr is deterministic, and thus one result is presented for benchmarks Navl 
and Nav2 (NavO was not tested by mm-HSolver's authors [21] )■ The mm-HSolver 
results are those reported in the literature. SpaceEx was run in deterministic 
mode on NavO (specifically, we set parameter 'directions' = 'box' [2"2"]). 

Parameter setting: We set the test duration T = 12sec, which we esti- 
mate is long enough to produce a falsifying trajectory for NavO if one exists. 
For SA+RED, we chose to generate 10 samples (|@| = 10). We will see that 
even this small number is enough for the algorithm to be competitive. A re- 
gression consists of 20 jobs. The SpaceEx parameters were varied in such a way 
that the approximation R of the reachable set R became increasingly precise. 
Clustering% was given the values 0, 20 and 80 (the smaller the Clustering%, the 
better the approximation and the longer the runtime) . The ODE solver timcstcp 
5 was given the values 0.0008,0.02,0.041 seconds. These are, respectively, the 
minimum, median, and average values of 5 used by the variable step-size ODE 
solver used by SA+RED. The smaller d, the better the approximation and the 
longer the runtime. The following parameters were fixed: 'directions' = 'box', 
'Local time horizon' = lOsec, rel-err = abs-err = 1.0e-10. The NavO SpaceEx 
configuration files can be obtained by request from the authors. 

The performance metrics: Each run produces a minimum robustness. 
For a given regression, we measure: the smallest, the average, and the largest 
minimum robustness found by the regression (min, avg, max in Table [1]). The 
standard deviation of minimum robustness is also reported (ct/). For SpaceEx, 
we had to simply assess whether R(j) intersected U or not. 

The cost metric: Each run also counts the number of simulated trajecto- 
ries in the course of its operation: SA simulates a trajectory for each proposed 
sample, SA+RED simulates a trajectory each time the constraint function of 
Prot(3[Wi] is evaluated (and for each sample), and mm-HSolver simulates tra- 



jectories in falsification mode. The trajectories simulated by SA and SA+RED 
have a common, fixed, prc-detcrmincd duration T . Thus the cost of these algo- 
rithms can be compared by looking at the Number of Trajectories (NT) each 
simulates (column NT in TableQ]- the overline denotes an average). The trajec- 
tories computed by mm-HSolver have varying lengths, determined by a quality 
estimate. So for comparison, we report the number of single simulation steps 
{SS), i.e. the number of points on a given trajectory (column SS - mm-HSolver, 
being deterministic, has one value of SS). Unfortunately, SS doesn't include the 
cost of doing verification in mm-HSolver, so it should be considered as a lower 
bound on its computational cost. On the other hand, because of the choice of 
T, the SS numbers reported for SA+RED should be treated as upper bounds: 
choosing a shorter a-priori T will naturally lead to smaller numbers. An exact 
comparison of the costs of SA+RED and mm-HSovler would require knowing 
the duration of the shortest falsifying trajectory, and setting the a-priori T to 
that, and somewhat incorporating the cost of verification. The operations that 
SpaceEx does are radically different from those of the other methods compared 
here. The only way to compare performance is through the runtime. 

Experimental setup: we impose an upper limit NTmax on NT: SA+RED 
is aborted when its NT reaches this maximum, and SA is made to generate 
NTmax samples. (Of course, SA+RED might converge before simulating all 
NTmax trajectories). 3 values were chosen for NTmax- 1000, 3000 and 5000. 
For each value, a regression is run and the results reported. This allows us to 
measure the competitiveness of the 2 algorithms (i.e. performance for cost). 

Experiments: Table Q] compares SA+RED to SA: we start by noting that 
SA+RED falsified Nav2, whereas SA failed to so. On most regressions, SA+RED 
achieves better performance metrics than SA, for the same (or lower) compu- 
tational cost. This is consistent whether considering best case (min), average 
case (avg) or worst case (max). There are 2 exceptions: for Navl and Nav2, 
NTmax = 5000 produces better average and max results for SA than for 
SA+RED. When running realistic system models, trajectory simulation is the 
biggest time consumer, so effectively NT is the limiting factor. So we argue that 
these 2 exceptions don't invalidate the superiority of SA+RED as they occur for 
high values of NT that might not be practical with real-world models. (In these 
cases, we observed that SA eventually produces a sequence of samples whose 
trajectories finish by making a large number of jumps between locations 3 and 
2, with a relatively high robustness. From there SA then produces a sample 
with (or close to 0) robustness. This happens on every Navl run we tried, and 
most Nav2 runs, resulting in the numbers reported. The RED step in SA+RED 
seems to avoid these trajectories by 'escaping' into local minima, and is worthy 
of further study.) 

Table H compares SA+RED to mm-HSolver. We note that SA+RED falsi- 
fies the benchmarks, as does mm-HSolver. For Navl, SS is greater than mm- 
HSolver's SS, though the falsifying runs have SS values (last column) both 
smaller and larger than mm-HSolver. For Nav2, which appears to be more chal- 



System 


/V 1 A/I A V 


NT 

((Jjvt) 


CT e 

J 


SA+RED Rob. 

min, avg, max 


SA Rob. 

min, avg, max 


NavO 


1000 
3000 
5000 


1004 (1.4) 
2716 (651) 
4220 (802) 


0.022 
0.019 
0.009 


0.2852, 0.30,0.35 
0.2852,0.29,0.32 
0.285,0.28,0.32 


0.2853,0.33,0.33 
0.2858,0.31,0.36 
0.286,0.32,0.35 


Navl 


1000 
3000 
5000 


662 (399) 
1129 (1033) 
1723 (1770) 


0.21 
0.23 
0.23 


0,0.43,0.65 
0,0.39,0.65 
0,0.38,0.68 


0,0.96,1.88 
0,0.99,1.80 
0,0,0 


Nav2 


1000 
3000 
5000 


902 (246) 
1720 (1032) 
1726 (1482) 


0.32 
0.3 
0.27 


0,0.54,0.78 
0,0.53,0.83 
0,0.62,0.79 


0.3089,1.11,1.90 
0.3305,1.29,1.95 
0,0.002,0.01 


Fosc 


1000 
3000 
5000 


1000 (9.3) 
3000 (8.7) 
5000 (11) 


0.024 
0.024 
0.028 


0.162,0.206,0.251 
0.163,0.203,0.270 
0.167,0.193,0.258 


0.1666,0.216,0.271 
0.173,0.212,0.254 
0.185, 0.218, 0.245 



Table 1. Comparison of SA and SA+RED. To avoid clutter, Robustness values are 
reported to the first differing decimal, with a minimum of 2 decimals, a f is standard 
deviation of robustness for SA+RED. 



System 


NT MAX 


SS 

{(TSS) 


°S 


SA+RED Rob 
min, avg, max 


mm-HSolver 
Rob, NT, SS 


SS at min Rob 
for SA+RED 


Navl 


1000 
3000 
5000 


47k (30k) 
79k (76k) 
143k (141k) 


0.21 
0.23 
0.23 


0,0.43,0.65 
0,0.39,0.65 
0,0.38,0.68 


0, 22,5454 


0,1560,16k 
0,0,1600, 127k 
7660, 38k, 102k, 159k 


Nav2 


1000 
3000 
5000 


63k (18k) 
126k (80k) 
124k (114k) 


0.32 
0.3 
0.27 


0,0.54,0.78 
0,0.53,0.83 
0,0.622,0.79 


0, 506, 138k 


2888, 74k 
14k, 57k, 210k 
3450, 121k, 331k 



Table 2. Comparison of SA+RED and mm-HSolver. The last column shows some of 
the SS values at which min robustness is achieved by SA+RED on various runs. 



lenging, SA+RED performed better on average than mm-HSolver. However, we 
point out again that exact comparison is hard. 

For SpaccEx running on NavO, we observed that our initial parameter set 
produces an R(j) that intersects U. Since this is inconclusive, we modified the 
parameters to ge t a be tter approximation. For parameter values (Clustering%, 
5) = (0,0.0008), R(j) and U were almost tangent, but SpaceEx runtimes far 
exceeded those of SA+RED (more than 1.5 hours). Moreover, SpaceEx did not 
reach a fixed point of its iterations (we tried up to j = 200 iterations before 
stopping due to high runtimes). Thus, we can not be sure that all of the reachable 
space was covered. While this may be seen as an analogous problem to the choice 
of T in SA+RED, the computational cost of increasing j is much more prohibitive 
than that of increasing T. We now present some detailed runtime results. For 
SA+RED, 'runtime' means the User time reported by the Unix time utility. 
SA+RED was run on a dedicated Intel Xeon processor, x86-64 architecture, 
under the Unix OS. SpaceEx reports its own runtime. It was run on a Dual- 



Clustering% 


8 (sec) 




SA+RED Runtime (sec) 


NTmax 




0.0008 


0.002 


0.041 


min,avg,max 




80 


737 


30 


15 


324, 426, 596 


1000 


20 


1066 


53 


33 


620, 1132, 1385 


3000 


10 


1460 


NA 


NA 


767,1617, 2216 


5000 





> 5400 


NA 


NA 







Table 3. Comparison of SA+RED and SpaceEx runtimes. NA means the experiment 
was not run, because a more accurate run was required. The right-most columns shows 
the NTmax constraint for which the SA+RED runtimes were obtained. 



Core Intel Centrino processor, under a Windows7 64b OS, with no other user 
applications running. 

Thus we may conclude that stochastic falsification and reachability analysis 
can play complementary roles in good design practice: first, stochastic falsifi- 
cation computes the robustness of the system with respect to some unsafe set. 
Guided by this, the designer may make the system more robust, which effectively 
increases the distance between the (unknown) reachable set and the unsafe set. 
Then the designer can run a reachability analysis algorithm where coarse over- 
approximations can yield conclusive results. 

5 Conclusions 

The minimum robustness of a hybrid system is an important indicator of how safe 
it is. In this paper, we presented an algorithm for computing a local minimum 
of the robustness for a certain class of linear hybrid systems. The algorithm can 
also be used to minimize the robustness of non-hybrid linear dynamic systems. 
When integrated with a higher-level stochastic search algorithm, the proposed 
algorithm has been shown to perform better than existing methods on literature 
benchmarks, and to complement reachability analysis. We will next deploy this 
capability to perform local descent search for the falsification of arbitrary linear 
temporal logic specifications, not only safety specifications. This investigation 
opens the way to several interesting research questions. Most practically, reduc- 
ing the number of tests NT results in an immediate reduction of the computation 
cost. Also useful, is the determination of an appropriate test duration T, rather 
than a fixed arbitrary value. 

In terms of performance guarantees, obtaining a lower bound on the optimum 
achieved in Problem [3] could lead to a lower bound on the optimal robustness. 
One level higher in the algorithm, it is important to get a theoretical under- 
standing of the behavior of the Markov chains iterated by SA+RED to further 
improve it. 
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